// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//   The Division Algorithm in a Multivariate Polynomial
//   Ring over a field
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//  This is the algorithm of IVA  2.3 Theorem 3.

//  First the Polynomial ring and some accessories.

Fld := GF(2);
n := 2;			  // number of variables
P<x,y> := PolynomialRing(F,n, "grevlex");
	// When not explicitly defined, Magma uses the lex term ordering,
	// with P.i >P.(i+1), P.i being the ith variable in P.

//  A home for the exponent of a monomial.
E := RSpace(Integers(), n);

// Some shorthand for Magma functions, relative to the term ordering,
// and a function for computing the leading exponent in the space E.

LT := LeadingTerm;
LC := LeadingCoefficient;
LM := LeadingMonomial;
LE := func < f | E ! [ Degree(LT(f),1), Degree(LT(f),2)] >;

//  A basic step in the division theorem is to check whether
//  a monomial divides another.  
//  Here are two versions of a function to check this.
//  The second uses some slick Magma tools.

DoesDivide := function(m, mm)
a := LE(mm) - LE(m);
num := Degree(a);
for i in [1..num] do
  if a[i] lt 0 then return false; end if;
end for;
return true;
end function;

DoesDivide := func < m, mm | forall{i : i in [1..n] | a[i] gt 0}
                           where a is LE(mm)- LE(m) >;
		// A slicker version of the same function.

// -----------------------------------------------
// Division Algorithm


// Input: sequence of polynomials f_i.
F  := [x^2 + x,  y^2 + y];
// Input: the polynomial we reduce modulo the f_i.
g := x^3* y^2;

// Output:  The quotients a_i and remainder r so that 
//    f = sum a_i f_i + r, and no terms of r are divisible 
//    by any of the LT(f_i). 

//    Here we initialize to 0.
A :=  [P ! 0 : i in [1..#F] ];
r  := P ! 0;

// A local copy of g that we be will altered until it is 0.
p := g;
while p ne  P!0 do
  i := 1;
  divisionocurred := false;
  while i le #F and divisionocurred eq false do	
     if DoesDivide(LM(F[i]), LM(p)) then
	A[i] := A[i] + (LT(p) div LT(F[i]));
	p     :=  p  - (LT(p) div LT(F[i])) * F[i];
        divisionocurred := true;
     else
	i := i+1;
     end if;
  end while;
  if divisionocurred eq false then
     r := r + LT(p);
     p := p - LT(p);
  end if;
end while;


//  Now you can check the result.
&+ [ A[i] * F[i]  : i in [1..#F] ] + r;
g;

//  Exercise:  Rewrite the algorithm as a function.

//  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//    Other Monomial Orderings
//  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Alternate orderings are "lex", "glex", "grevlex", 
// "elim", "univ", "weight".
// "elim" may be used in two ways.  
//  In the first, it takes a parameter k in [1..n].
//  In the second it takes a parameter U a subsequence of [1..n].
//  The first is equivalent to U= [1..k].
// "univ" takes a parameter k.
//  It is equivalent to "elim" with U = Exclude([1,..n], k);
//  You can also use variable names instead of the index.
// Examples:
k := 2;               // You can use k in [1..6] here
P<s,t,x,y,z,w> := PolynomialRing(F, 6, "elim", k);

U := [1,3,5];         // Any subset U of [1..6] or [s,t,x,y,z,w]
P<s,t,x,y,z,w> := PolynomialRing(F, 6, "elim", U);

P<s,t,x,y,z,w> := PolynomialRing(F, 6, "elim", k);

// "weight" takes a parameter of n^2 nonnegative rational numbers.
//  These are parsed by Magma into n vectors W_i, i in [1..6].
//  The monomial order is that given by the W_i
//  This gives graded lex with r < s< x < y < z < w
W := [ 1,1,1,1,1,1, 0,1,1,1,1,1, 0,0,1,1,1,1, 
       0,0,0,1,1,1, 0,0,0,0,1,1, 0,0,0,0,0,1 ];
P<s,t,x,y,z,w> := PolynomialRing(F, 6, "weight", W);