-- -*- coding: utf-8 -*-
-- In this tutorial we introduce a number of basic operations 
-- using Gröbner bases, and at the same time become familiar
-- with a range of useful Macaulay2 constructs.

----------------------------
-- A. A first Gröbner basis 
----------------------------

-- To compute the Gröbner basis of an ideal
-- $(x^2y,xy^2+x^3)$ in the polynomial ring in
-- four variables we proceed as follows:
--
-- Our favorite field
KK = ZZ/32003

-- The polynomial ring

R = KK[x,y,z,w]

-- and the ideal

I = ideal(x^2*y,x*y^2+x^3)

-- now the punch line:

J = gens gb I

-- This is the Gröbner basis of $I$ under the graded
-- reverse lexicographic order.  In Macaulay2, monomial
-- orders are associated with a polynomial ring.

-- For example, the lexicographic order is specified using:
R = KK[x,y,z,w,MonomialOrder=>Lex]
I = substitute(I,R)
gens gb I

-- The Gröbner basis is the same, since this is a small 
-- example.  The polynomials are sorted in ascending monomial order
-- by their lead terms, but otherwise the two Gröbner bases are
-- the same here.

------------------------------------------------
-- B. Random regular sequences
------------------------------------------------

-- An interesting and illustrative open problem
-- is to understand the initial ideal (and 
-- the Gröbner basis) of a ``generic'' 
-- regular sequence.  To study a very simple case
-- we take a matrix of 2 random forms 
-- in a polynomial ring in
-- 3 variables:

R = KK[x,y,z]
F = random(R^1, R^{-2,-3})

-- makes $F$ into a $1 \times 2$ matrix whose elements
-- have degrees $2,3$ (that is, $F$ is a random map
-- to the free module $R^1$, which has its one
-- generator in the (default) degree, $0$, from
-- the free module with generators in the listed
-- degrees, $\{2,3\}$).  We now can compute

GB = gens gb F
LT = leadTerm GB
betti LT
-- betti is a routine that displays degrees of generators
-- and also in free resolutions (which we will learn about later).

-- In this case, the output
-- shows that there are Gröbner basis elements
-- of degrees 2,3, and 4.  This result is
-- dependent on the monomial order in the ring $R$;
-- for example we could take the lexicographic
-- order

R = KK[x,y,z, MonomialOrder => Lex]

-- (see {\tt help MonomialOrder} for other possibilities).

-- We get 
F = random(R^1, R^{-2,-3})
GB = gens gb F
LT = leadTerm GB
betti LT

-- and there are Gröbner basis elements of degrees 
-- $2,3,4,5,6.$

-----------------------------------------------
-- C. Division With Remainder
-----------------------------------------------
-- A major application of Gröbner bases is
-- to decide whether an element is in a given
-- ideal, and whether two elements reduce to
-- the same thing modulo an ideal.  For
-- example, everyone knows that the trace
-- of a nilpotent matrix is 0. We can produce
-- an ideal $I$ that defines the variety $X$ of 
-- nilpotent $3 \times 3$ matrices by taking the cube
-- of a generic matrix and setting the entries
-- equal to zero.  Here's how:

R = KK[a..i]
M = genericMatrix(R,a,3,3)
N = M^3
I = flatten N

-- (actually this produces a 1 x 9 matrix of
-- of forms, not the ideal: {\tt J = ideal I};
-- the matrix will be more useful to us).
-- But the trace is not in $I$!  This is obvious
-- from the fact that the trace has degree $1$,
-- but the polynomials in $I$ are of degree $3$.
-- We could also check by division with
-- remainder:

Tr = trace M 
Tr //I  -- the quotient, which is 0
Tr % I  -- the remainder, which is Tr again

-- (Here {\tt Tr} is an element of $R$, not a matrix.
-- We could do the same thing with a $1 \times 1$ matrix
-- with {\tt Tr} as its element.)
-- This is of course because the entries of $I$ do
-- NOT
-- generate the ideal of all forms 
-- vanishing on $X$ -- this we may find with
-- {\tt J = radical ideal I},
-- (but this takes a while: see the documentation for
-- {\tt radical} on a faster way to find this)
-- which shows that the radical is generated by
-- the trace, the determinant, and the sum of 
-- the principal $2 \times 2$ minors, that is, by the
-- coefficients of the characteristic polynomial.
-- In particular, we can try the powers of the
-- radical:

Tr^2 % I
Tr^3 % I
Tr^4 % I
Tr^5 % I
Tr^6 % I
Tr^7 % I

-- The seventh power is the first one in the 
-- ideal!  (Bernard Mourrain has worked out a
-- formula for which power in general.)
-- In this case

Tr^6 // I

-- is not 0.  It is a matrix that makes the
-- following true:

Tr^6 == I * (Tr^6 // I) + (Tr^6 % I)

----------------------------------------------
-- D. Elimination Theory
----------------------------------------------

-- Computing the elimination ideal $I \cap k[xi, \ldots, xn]$
-- is one of the most important applications of Gröbner bases.

R = KK[t,y,z,MonomialOrder=>Lex]
I = ideal(y-(t^2+t+1), z-(t^3+1))
GB = gens gb I
F = GB_(0,0)
substitute(F, {y =>t^2+t+1, z=>t^3+1})

-- F is the polynomial that gives an algebraic relation between
-- $t^2+t+1$ and $t^3+1$.

-- Another way to accomplish this in Macaulay2 is to use the
-- {\tt eliminate} function.  In this case, the monomial order of
-- the ring is not important.
R = KK[y,z,t]
I = substitute(I,R)
eliminate(I,t)

-- Yet another method is to use ring maps.
A = KK[t]
B = KK[y,z]
G = map(A,B,{t^2+t+1, t^3+1})
kernel G

---------------------------------------
-- E. Intersections and ideal quotients 
---------------------------------------

-- An interesting problem is to investigate ideals of the
-- form $I = (x^d, y^d, z^d, f)$, where $f$ is a polynomial,
-- in three variables $x,y,z$.
-- 
-- First, let's compute the ideal quotient for an example,
-- by using the method given in class.
-- 
R = KK[t,x,y,z]
I = ideal(x^3,y^3,z^3)
F = x+y+z
L = t*I + (1-t)*ideal(F)
L1 = eliminate(L,t)
gens gb L1
-- Each of these is divisible by F:
(gens L1) % F
-- Divide by F.
J = ideal ((gens L1)//F)
mingens J
betti oo

-- J is defined by 5 equations: the original 3, another cubic
-- and also a quartic. (oo is the previous output).

-- Now, let's do this the easier way.
R = KK[x,y,z]
I = ideal(x^3,y^3,z^3)
F = x+y+z
J = I : F
betti J
transpose gens J
transpose gens gb J

-- Problem: how many generators can you obtain using this construction,
-- where you may try different positive integers $d$, and polynomials $f$?

----------------
-- F. Saturation
----------------

-- A very useful construction is the saturation of an ideal
-- $I$ with respect to a polynomial $f$, or another ideal $J$.
R = KK[t,a..f]
I = ideal(a*b*c-d*e*f, a^2*b-c^2*d, a*f^2-d*b*c)
F = a*b*c*d*e*f
J = eliminate(I + ideal(t*F-1), t)
transpose gens J

-- There is a builtin routine in Macaulay2 for computing saturations:
R = KK[a..f]
I = substitute(I,R)
F = product gens R
J' = saturate(I,F)
transpose gens J'