##############################################################################
### variable ordering
##############################################################################
`varorder/get` := proc()
	global  `varorder/ordering`, `varorder/vars`; 
  option inline;
  `varorder/ordering`(`varorder/vars`);
end:

varorder := proc()
  global `varorder/vars`, `varorder/ordering`, `type/rtrm`, `type/rmon`, `type/rpol`, `type/rcc`;
      `type/rtrm` := monomial(rational, {`varorder/vars`}):
      `type/rmon` := monomial(identical(1), {`varorder/vars`}):
      `type/rpol` := polynom(anything, {`varorder/vars`}): 
      `type/rcc`  := [rpol, rmon, set] :# set(`varorder/vars`)
  if nargs = 0 then 
    return (`varorder/get`());
  elif nargs = 2 then  
    if type(args[1], symbol) and type(args[2], list(symbol)) then
      `varorder/ordering` := args[1];
      `varorder/vars` := op(args[2]);
      procname();
    else
      error "Wrong type(s) of arguments";
    fi;  
  else
    error "Wrong number of arguments";
  fi;  
end:

##############################################################################
### polynomial utilities
##############################################################################

if StringTools[IsPrefix]("Maple 9.", kernelopts(version)) then  # Maple 9.x

  LMon := proc(p::rpol) option inline; Groebner[leadterm](p,`varorder/get`()); end:
  LCoe := proc(p::rpol) option inline; Groebner[leadcoeff](p,  `varorder/get`()); end:
  LTrm := proc(p::rpol) option inline; `*`(Groebner[leadmon](p, `varorder/get`())) end:
  RPol := proc(p::rpol) option inline; p - LTrm(p) end:
  
  CoeMon := proc(m::rtrm)  
    # decomposes term c*m to (c, m)
    option inline;
    Groebner[leadmon](m, `varorder/get`()) ### TODO: asi neefektivni
  end:
  
  # sort
  RSrt := proc(T)
    #option inline; ### cannot be inlined!
    sort(T,(t1,t2)->Groebner[testorder](t2,t1,`varorder/get`())) # reversed order
  end:

elif assigned(Groebner[LeadingMonomial]) then  # Maple >=10

  LMon := proc(p::rpol) option inline; Groebner[LeadingMonomial](p,`varorder/get`()); end:
  LCoe := proc(p::rpol) option inline; Groebner[LeadingCoefficient](p,  `varorder/get`()); end:
  LTrm := proc(p::rpol) option inline; `*`(Groebner[LeadingTerm](p, `varorder/get`())) end:
  RPol := proc(p::rpol) option inline; p - LTrm(p) end:
  
  CoeMon := proc(m::rtrm)  
    # decomposes term c*m to (c, m)
    option inline;
    Groebner[LeadingTerm](m, `varorder/get`()) ### TODO: asi neefektivni
  end:

  # sort
  RSrt := proc(T)
    #option inline; ### cannot be inlined!
    sort(T,(t1,t2)->Groebner[TestOrder](t2,t1,`varorder/get`())) # reversed order
  end:
fi:

`rpol/decomp` := proc(p::rpol)
  # decompose polynom to a sequence of its monoms
  option inline;
  `if`(type(p,rtrm), p, op(p)); 
end:

# collect
RCol := proc(p::rpol)::rpol; 
	global  `varorder/vars`; 
  option inline;
  #expand(p);
  collect(p, [`varorder/vars`], distributed) 
end:

# extremal term (returns the term with mononom of extreme ranking)
RMax := proc(s) option inline; RSrt(convert(s,list))[1] end: ### TODO: neefektivni
RMin := proc(s) option inline; RSrt(convert(s,list))[nops(s)] end: ### TODO: neefektivni

# decompose polynom to substitution
spol := proc(p)
    local sp, lm, rp, lc, rp1;
    if type(p,equation) then 
	    sp := RCol(lhs(p) - rhs(p) )
	  else
	    sp := RCol(p);
	  fi;
	  if sp=0 then
      return NULL;
    else
	    lm :=  LMon(sp);
	    rp :=  RPol(sp); 
	    lc :=  LCoe(sp);
	    rp1 := -rp/lc;
	    return (lm=rp1);
	  fi;
end:

# number of terms in polynom
mlen := proc(p::rpol)
  option inline;
  `if`(type(p,rtrm), 1, nops(remove(type, p, constant)));
end:

##############################################################################
### linear algebra
##############################################################################

CoeffMatrix := proc(Ps::sequential(rpol))
  # generalization of LinearAlgebra[GenerateMatrix]
  # set of all monoms in a set of polynomials Ps are considered as variables
  # we collect the coefficients
  local C, ms, pn, mn, i, t, c, m, Ms;
  ms := [AllMonoms(Ps)];
  mn := nops(ms);  
  Ms := table([seq(ms[i]=i, i=1..mn)]); # [monom = its position in ms, ...]    
  pn := nops(Ps);
  C := Matrix(storage=sparse, pn, mn);
  for i from 1 to pn do
    # store coefficients of Ps[i] to the table
    for t in [`rpol/decomp`(Ps[i])] do
      c, m := CoeMon(t);
      #print(m, [i,Ms[m]] = c);
      C[i,Ms[m]] := c;
    od;
  od; 
  return(C, ms ); # C = matrix of coefficients, ms = vector of variables
end:

AllMonoms := proc(p)
  option inline;
  op(RSrt([op(`AllMonoms/s`(args))]))
end:

`AllMonoms/s` := proc(p) 
  if nargs > 1 then 
    `union`(op(map(procname, [args])))
  else
    if type(p, sequential) then
      `union`(op(map(procname, p)))
    elif type(p, rpol) then
      `AllMonoms/1`(p)
    else
      error "Wrong type of argument, must be rpol", p;
    fi
  fi
end:

`AllMonoms/1` := proc(p::rpol)
  option inline;
  map(LMon, {`rpol/decomp`(p)});
end:

##############################################################################
### reducetable
##############################################################################

`type/reducetable`:= table:

`reducetable/indices` := proc (R::reducetable)
   option inline;
   map(op, [indices(R)]) 
end :

`reducetable/getsubs` := proc(R::reducetable, rs::sequential(rmon))
  if nops(rs) > 0 then 
      return op(Map(m -> m =  R[m], rs));
  fi;
end:

`reducetable/add` := proc(R::reducetable, Rm::uneval, s::{rpol, sequential(rpol)})
	local sp, lm, rp, lc, rp1, cc;
	if nargs>3 then
		return op(map(s1->`reducetable/add`(R,Rm,s), [args][3..-1]));
	elif type(s, sequential) then
		return map(s1->`reducetable/add`(R, Rm, s1), s);
	elif type(s,equation) then 
	  sp := RCol(lhs(s) - rhs(s) )
	elif type(s,rpol) then
	  sp := s; # RCol(s);
	else
		error "Unsupported type of parameter", s;
	fi;

	lm :=  LMon(sp);
	if type(sp,constant) or lm=0 then 
	  error "Trying to add a constant to the reducetable", s;
	fi;
	rp :=  RPol(sp); 
	lc :=  LCoe(sp);
	rp1 := -rp/lc;
	
	if assigned(R[lm]) then 
	  printf( "`reducetable/add`:\n");
	  printf( "During adding s-polynom %q to the reducetable\n%q\n", sp);
	  printf( "Trying to assign to existing item of index %q",lm); 
	  printf( " new value %q, ", normal(rp1));
	  printf( " old value %q.\n", normal(R[lm]));
	  error "Trying to assign to existing item of index ",lm,"to the reducetable";  
	fi;
	
	userinfo(5, reduce, 
	  "Adding ", sp ,  "to reducetable into item ",lm, " value", rp1);
	
  #if mlen(rp1) <= 1 then
  #   `reducetable/reduce`(R, lm=rp1); ### ?
  #fi;
  
  forget(`reduce`);
  forget(`reduce/1`);
  forget(`cc/single2`);
  
  push(Rm, lm);
  assign(R[lm] , rp1); ### operands to assign 
   
end:

`reducetable/transform` := proc(R, Rm::uneval)
	local sps, M, ms, N, R1, rs;
	if type(R, reducetable) then 
		sps := `reducetable/polynomlist`(R); # get a list of all S-polynomials
	elif type(R, sequential(rpol)) then
	  sps := convert(R, list);
	else
		error "Unsupported parameter type";
	fi;
	M, ms := CoeffMatrix(sps);
	N  := LinearAlgebra[GaussianElimination](M);
	R1 := table();
	rs := convert(evalm(N &* ms), list);
	#rs :=  remove(type, rs, constant);
	`reducetable/add`(R1, Rm, rs);
	return(eval(R1));
end:

`reducetable/polynomlist` := proc(R::reducetable)
	 map(sp -> lhs(sp) - rhs(sp), op(op(R)));
end:

`reducetable/reduce` := proc(R::reducetable, s)
  # reduce all items in the table by given S-polynom
  local rsp;
  map(proc(sp)
        R[lhs(sp)] :=  Asubs(s, rhs(sp));
      end, op(op(R)));
  op(R);
end:


##############################################################################
### reduce 
##############################################################################

reduce := proc (p::rpol, R::reducetable, Rm::uneval, RCache::table)
  # reduce polynom by table 
  option remember, system;  
  local p0, p1, rp, lm, lc, lmr;
  if assigned(RCache[p]) then 
  	p0 := RCache[p];
  else
    p0 := p;
  fi;
  RCache[p] := `reduce/1`(p0,R, Rm);
end:

`reduce/1` := proc(p::rpol, R::reducetable, Rm::uneval)
  local p1, rp, lm, lc, lmr;
  p1 := p;#RCol(p0);
  reportf(5, "`reduce/p`: reducing %a=%a\n", p, p1); 
  if type(p1, constant) then 
    return p1;
  else
    lm := LMon(p1);
    lc := LCoe(p1);
    rp := RPol(p1);
    reportl(4, "`reduce/p`: reducing", (p = p1) = ({lc} * {lm}) + rp);
    lmr := `reduce/mon`(lm, R, Rm);
    if lmr <> FAIL then
      procname(RCol(lc*lmr) + rp, R, Rm);
    else
      RCol(lc*lm) + procname(rp, R, Rm);
    fi; 
  fi;
end:

`reduce/mon` := proc(m::rmon, R::reducetable, Rm)
  local r,q;
  r, q := `reduce/getsubs`(m, R, Rm);
  if r <> FAIL then  q*R[r] else FAIL fi;
end:

`reduce/getsubs` := proc(m::rmon, R::reducetable, Rm)
  option remember; # mandatory
  local r, q;
  #for r in eval(Rm) do
  for r in `reducetable/indices`(R) do
   if divide(m,r,q) then
       return (r, q);
   fi  
  od:
  FAIL, FAIL;
end: 


##############################################################################
### gbasis 
##############################################################################


gbasis := proc(F)
  local Rm, R, c, c1, rc, erc, erc1, erc2, p,  m, M, ms, N,H, s, sizedc, konec, RCache;
  H := table();
  Rm := [];
  R := `reducetable/transform`(F, Rm);
  RCache := table();
  reportf(2, "gbasis: Starting with %q\n", op(R), Rm);
  konec := false;
  while konec=false do
    c := cc(Rm, R, H);
    if nops(c)=0 then
      reportf(2, "\ngbasis: END, no cc.\n");
    	konec := true;
    else
      c := `cc/sort`(c);
      s := 1;
      rc := NULL;
      reportf(2, "\ngbasis: Selecting from cc %q\n", c);
      while nops(c)>0 and s < 15 do 
      	p := 0;
      	while nops(c)>0 and p=0 do
          c1 := pop(c);
          if assigned(H[c1[2]]) then 
            H[c1[2]] := [H[c1[2]][1], H[c1[2]][2] union c1[3]];
          fi;
          # reduce s-polynomial from c1
          p := reduce(c1[1], R, Rm, RCache);
        od;  
        if p<>0 then
        	s := s*mlen(p);
          rc := rc, p;
        fi;
      od;	      
      
      erc := [];  
      if rc = NULL and nops(c)=0 then
      	reportf(1, "gbasis: END, all cc reduced to 0.\n");
      	konec := true;
      elif nops([rc])=1 then
        erc := [rc];
        reportf(2, "gbasis: Selected cc is %q\n", erc);
      else  
      	reportf(2, "gbasis: Selected cc are %q\n", rc);
        M, ms := CoeffMatrix([rc]);
        N := LinearAlgebra[GaussianElimination](M);
        erc := convert(evalm(N &* ms), list);
        erc := remove(is, erc, 0);  
        reportf(1, "gbasis: Eliminated cc are %q\n", erc);	
      fi;
      
      if nops(erc)>0 then
        erc1, erc2 := selectremove(type, erc, constant);
        if nops(erc1)>0 then
        	reportf(1, "gbasis: END, 1=0.\n");
          konec := 1=0;
        else	
         reportf(3, "gbasis: Adding cc %q\n", erc2);
         `reducetable/add`(R, Rm, erc2);
        fi;
      else
        reportf(1, "gbasis: END, all cc eliminated.\n");	
      	konec := true;
      fi;
  fi
  od;
  if konec = (1=0) then 
    {1}
  else
    {`reducetable/getsubs`(R,`monom/min`(`reducetable/indices`(R)))};
  fi;
end:


##############################################################################
### cc 
##############################################################################

cc := proc (M::sequential(rmon), R::reducetable, H::table)
	 #option inline;
	 reportf(4, "\n`cc`: computing cc of %q\n", M);
   map(`cc/m`, [cp(H, op(M))], R);
end:

`cc/m` := proc (m::{rmon=set(set)}, R::reducetable)
  local r, S1, SS, res;     
  r := lhs(m);
  S1 := rhs(m) [1];
  SS := rhs(m) [2..-1];

  res := op(map(S2 -> `cc/m/strategy`(R, r, S1, S2), SS));
  reportf(4, "`cc/m`: cc of %a is %q\n", m, res);

  return res;
end:

`cc/m/strategy` := `cc/m/all`:

`cc/m/single` := proc (R::reducetable, r::rmon, S1, S2)
    local res;
    res := `cc/sort`([`cc/m/all`(R, r, S1,S2)])[1];
    reportf(3, "`cc/m/single`: choosen cc=%q\n", res);
    res;  
end:

`cc/m/all` := proc (R::reducetable, r::rmon, S1, S2)
	  option inline;
	  op(mapmap((s1,s2) -> `cc/1`(R, r, s1, s2),
	             S1,S2)); 
end:	


`cc/1` := proc (R::reducetable, r::rmon, r1::rmon, r2::rmon)
  local  R1, R2, rs1, rs2;  
  #if r1 <> 1 then  rs1 := R[r1] else rs1 := 1 fi; 
  #if r2 <> 1 then  rs2 := R[r2] else rs2 := 1 fi;
  rs1 := R[r1];
  rs2 := R[r2];
  #reportf(2, "`cc/1`: cc of %a ---{%a, %a}---> {%a, %a} = \n", r, r/r1, r/r2,   r1, r2); 
  return [RCol(r/r1*rs1 - r/r2*rs2), r, indets(r/r1) union indets(r/r2)];  ### ? 
  # [s-pol, source monom, {source directions}]
end:

`cc/sort` := proc(cs::sequential(rcc))
  # sort cc by size
  local c, sizedc;
  sizedc := map(p -> p=mlen(p[1]), cs);  # add sizes
  c := sort(sizedc, (p,q)->rhs(p)<rhs(q)); # sort by sizes
  c := map(lhs,c); # forget sizes
  return c;
end:



##############################################################################
### MM 
##############################################################################

`type/count` := monomial(identical(1)):

`monom/divides` := proc(t::count, s::count)
  local q;
  evalb(divide(s,t,q) and q <> 1)
end:

`monom/min` := proc(S::{set(count),list(count)})
  remove(
    proc(s,S)
      local t;
      for t in S do if `monom/divides`(t,s) then RETURN(true) fi; false od
    end, S, S);
end:

`monom/hull` := proc()
  local M,N,A,MA; 
  global NT;
  M := {args};
  N := {};
  while M <> N do
    N := M;
    for A in M do
      MA := select(proc(B,A) evalb(A intersect B <> {}) end, M, A);
      M := M minus MA union {`union`(op(MA))};
    od
  od; 
  M
end:

`monom/sets` := proc(p::count, S::{set(count),list(count)})
  local s,b;
  b := [seq(p/s, s = S)];
  op(map(indets, select(type,b,count)));
end:

cp := proc(H::table)
  local M, aa, c1, c2, T1, T2, ms;
  aa := {args[2..-1]};
  M := `monom/min`(aa);
  c1 := `cp/I`(aa, M);;
  c2 := `cp/II`(H, aa, M);
  T1 := table(c1);
  T2 := table(c2);
  ms := map(lhs, c1) union map(lhs, c2);
  op(map(m ->m=`cp/u`(T1, T2,m) ,ms))
end:

`cp/u` := proc (T1,T2,m)
	option inline;
	`if`(assigned(T1[m]), 
	     `if`(assigned(T2[m]), 
	          T1[m] union T2[m], 
	          T1[m]), 
	      T2[m]);
end:

`cp/I` := proc(L::set, M::set)
  map(
    proc(u, M)
      u = {{u}, select(proc(v,u) `monom/divides`(v,u) end, M, u)};
    end, L minus M, M);
end:

`cp/II` := proc(H::table, L::set, M::set)
  local u,v,p,S,R;
  S := {seq(seq(proc(u,v) if u <> v then lcm(u,v) fi end(u,v), u=M), v=M)};
  #R := {seq(p =  `monom/hull`(`monom/sets`(p, M)), p = S)};
  R := {seq(p = `cp/II/hull`(H, p, M), p = S)};
  R := remove(proc(e) evalb(nops(rhs(e)) = 1) end, R);
  map(proc(e, L) lhs(e) = map(`cp/II/1`, rhs(e), lhs(e), L) end, R, L)
end:

`cp/II/hull` := proc(H::table, p, M) # HB
  local h, old, oldh, oldM;
  if not(assigned(H[p])) then
    h := `monom/hull`(`monom/sets`(p, M))
  else
    old := H[p];
    oldM := old[1];
    oldh := old[2];
    h := `monom/hull`( `monom/sets`(p, M minus oldM), (op(oldh)));
  fi;
  H[p] := [M, h];
  return h;
end:

`cp/II/1` := proc(rs::set, u::count, L::set)
  option inline;
  select(proc(s,rs,u) evalb(`monom/divides`(s,u) and (indets(u/s) subset rs)) end, L, rs, u)
end:

##############################################################################
### auxiliary
##############################################################################

mapmap := proc(f, A, B) 
  local ff; 
  ff := rcurry(f, args[4..-1]); 
  map(a-> op(map(b -> ff(a,b), B)), A) 
end:

push := proc(L::uneval, x) L := [op(eval(L)),x] end:
pop  := proc(L::uneval) local x; x := L[-1]; L := L[1..-2];  x; end:


Asubs := proc() inline; `if`(nargs<=1, args, aalgsubs(args)) end: 


aalgsubs := proc(p,f) inline; `if`(type(p, list), rasubs(args), algsubs(args)) end:

rasubs := proc(p,f)
  local n;
  n := nops(p);
  if n = 0 then f
  elif n=1 then algsubs(op(p), args[2..-1]) 
  else procname(p[2..-1], algsubs(p[1], args[2..-1]), args[3..-1]) fi;
end:  

report := proc(i::integer)
  if i=0 or (assigned(infolevel[reduce]) and infolevel[reduce] >= i) 
  	or (assigned(infolevel[all]) and infolevel[all] >= i)  then 
      print(args[2..-1]) 
  fi;
end:
reportl := proc(i::integer)
  if i=0 or  (assigned(infolevel[reduce]) and infolevel[reduce] >= i) 
  	or (assigned(infolevel[all]) and infolevel[all] >= i)  then 
      lprint(args[2..-1]) 
  fi;
end:
reportf := proc(i::integer)
  if i=0 or  (assigned(infolevel[reduce]) and infolevel[reduce] >= i) 
  	or (assigned(infolevel[all]) and infolevel[all] >= i)  then 
      printf(args[2..-1]) 
  fi;
end:

##############################################################################
### Maple backward compatibility
##############################################################################

# < Maple 11
if not assigned(Map) then Map := map fi: # paralell version of map 

# Maple V
if not assigned(`subset`) then 
  kernelopts(version);
  error "Unsupported Maple version.";
  `subset` := proc (A,B) option inline; evalb(A minus B = {}) end: 
fi;










